The relationships between microbiome diversity and epidemiology in domestic species of malaria-mediated mosquitoes of Korea

Microbiota in the mosquito plays an important role in their behavior and vector competence. The composition of their microbiome is strongly influenced by the environment, especially their habitat. The microbiome profiles of adult female Anopheles sinensis mosquitoes from malaria hyperendemic and hypoendemic areas in Republic of Korea were compared using 16S rRNA Illumina sequencing. In different epidemiology groups, the alpha and beta diversity analyses were significant. The major bacterial phylum was Proteobacteria. The most abundant species in the microbiome of hyperendemic mosquitoes were the genera Staphylococcus, Erwinia, Serratia, and Pantoea. Notably, a distinct microbiome profile characterized by the dominance of Pseudomonas synxantha was identified in the hypoendemic area, suggesting a potential correlation between the microbiome profiles and the incidence of malaria cases.

Data at the genus level were analyzed to identify bacterial operational taxonomic units (OTUs) present in all mosquito samples. (Fig. 3C,D, Supplementary Information 3). A total of 51 genera were found, the five most abundant bacterial genera (average abundance) were Pseudomonas (29.22%), Staphylococcus (10.5%), Erwinia (9.37%), Serratia (6.91%) and Acinetobacter (6.7%) (Fig. 3C,D, Supplementary Information 3). Variations in microbiome profiles were observed across different regions within the hyperendemic area. The hyperendemic area 1 nearest to the northern border had a greater Staphylococcus ratio than the other regions (Figs. 1, 3, Supplementary Information 3). Brevibacterium was found in high abundance in hyperendemic area 1S, however, this microbiome was not found in any other regions except in hyperendemic area 6S. Staphylococcus was also found at a higher rate than in other regions in hyperendemic area 3, which is the closest region to hyperendemic area 1. In hyperendemic area 3, Erwinia and Enterobacteriaceae_uc dominate the microbiome profile. Erwinia was also observed in a high ratio in the hyperendemic areas 4M and 7M. Despite being close-range regions, hyperendemic areas 4 and 5 showed distinct microbiome profiles, and there were also substantial differences between organs. Acinetobacter dominated the microbiome in hyperendemic area 4S, followed by Chryseobacterium, Erwinia, and Pseudomonas. The microbiome profile of the hyperendemic area 5M was dominated by Enterobacteriaceae_uc, followed by Serratia. Acinetobacter dominated the microbiome profile in hyperendemic area 5S, followed by Arcobacter and Pantoea. In hyperendemic area 6, the mid-gut and salivary glands had a microbiome profile dominated by Serratia and Gibbsiella. The microbiome profile of hyperendemic area 7 varied by organ. Pantoea and Erwinia were dominant in Hyperendemic area 7M, whereas Staphylococcus and Asaia were dominant in Hyperendemic area 7S. Hyperendemic area 8S showed a distinct microbiome profile dominated by Arcobacter. Hyperendemic areas 8M and 2M had similar microbiome to hypoendemic areas. Pseudomonas dominated microbiome profile from all hypoendemic areas ( Supplementary Information 3).

Sample
Taxon name p-value p-value (FDR) Taxonomic relative abundance (%)  There are 21 microbiotas shared by the mid-gut and salivary glands, as well as 20 and 10 organ-specific microbiotas, respectively. The microbiome shared by both organs had a high ratio, and Arcobacter was the only organ-specific microbiome in the salivary glands with a percentage of above 5%. In Acinetobacter, there was a statistically significant difference (Wilcoxon rank-sum test, p < 0.05). www.nature.com/scientificreports/ Epidemiological differences in mosquito microbiome composition. The data was further compared based on patient incidence, which allowed for a division between hyperendemic and hypoendemic areas.
There was a significant difference in alpha and beta diversity analyses based on epidemiology (Table 3, Figs.2, 6). The hyperendemic and hypoendemic areas were significant in PERMANOVA analysis using the Jaccard distance, which evaluates dissimilarity across data sets (Table 3). Furthermore, in both areas, the Shannon diversity index, which is a measure of the number (abundance) and relative abundance (uniformity) of species in the ecosystem, was significant (Fig. 2). Both areas shared 11 microbiotas mainly Pseudomonas and Acinetobacter. According to the analysis conducted, mosquitoes from hyperendemic areas exhibited 17 unique taxa, which were dominated by Staphylococcus (15.56%), Erwinia (14.05%), Pseudomonas (10.46%), Serratia (10.25%), Acinetobacter (9.48%) (Fig. 5). In hypoendemic areas, Pseudomonas (66.73%) dominates the microbiome profile, followed by Anaerobacillus (4.49%) (Fig. 5B). The microbiome profile of hyperendemic and hypoendemic areas have significant differences on Pseudomonas and Pantoea (Fig. 5B,C). Pantoea was only observed in hyperendemic areas (Fig. 5B,C). Pseudomonas was found at a low proportion in the hyperendemic area but was dominating in the hypoendemic area (Fig. 5B,C).
Principal coordinated analysis (PCoA) analysis showed great variance in microbial communities depending on epidemiology and collecting site (Fig. 6). In the graph, except for hyperendemic area 2-M, which is close to the sample of the hypoendemic area in graph, the hyperendemic area and hypoendemic area showed different patterns. Samples in the hypoendemic area were clustered regardless of region, but samples in the hyperendemic area had a clustering pattern according to region. In the hyperendemic area, 3 clustering were confirmed (Group

Discussion
The definition of microbiome per se has recently expanded from microorganisms interacting with a host to the total of microbiota and their gene products related to environment and hosts 7 influencing biology, longevity, behaviors, nutrition, and immunity [25][26][27] . The host's microbiota has the potential to co-evolve and continually affect the host's health 8 . There are several different environments where the microbiota may be acquired, including through food intake and wound invasion. As a result, the habitat of the host can affect the microbiome of the mosquito and may change based on regional environmental factors 27 . This study aimed to investigate whether there are significant differences in the microbiome profiles of mosquitoes between hyperendemic and hypoendemic areas of malaria incidence. The main goal of this study was to identify microbe species-biomarkers in different endemic areas which involved in malaria incidence. The current study demonstrates that the microbiome profile of mosquitoes collected in malaria-hyperendemic and hypoendemic areas was divided into distinct patterns at the level of genus and species. Mosquitoes in hyperendemic areas contained area-specific microbiome such as Serratia, Staphylococcus, Erwinia, Enterobacter, and Arcobacter, while Pseudomonas is dominant in hypoendemic areas, implying that malaria incidence can be analyzed with area-specific biomarkers. Remarkably, our findings  Table 3. An. sinensis microbiome pairwise beta diversity comparisons from different types of collecting sites. The comparison of Jaccard distance in hyperendemic and hypoendemic areas revealed significant differences. However, The comparison of Bray-Curtis distance in hyperendemic and hypoendemic areas has no significant differences. PERMANOVA (999 permutations) tests with Benjamini-Hochberg FDR correction were used to make these comparisons (q-value). The significance level was set at q < 0.05 , with n = the number of pools processed and each pool containing 5 individual mosquitoes. www.nature.com/scientificreports/ revealed that the heterogeneity of microbiome profiles was higher in the hyperendemic area compared to the hypoendemic area, despite the spatial distance between hyperendemic areas being closer than in the hypoendemic area. It is suggested that malaria incidence in ROK might be related to northern border areas of malaria hot zones such as North Korea where malaria cases were reported in over 5000 patients in 2016 28 . As a result of human migration, cross-border population transfers, ecological changes, and vector population dynamics, border malaria frequently occurs in some transmission zones 29 . Malaria outbreaks on the Thai-Myanmar and Thai-Cambodia borders are typical border outbreaks, and they are caused by constant population migration, the movement of malaria patients, abuse of self-medication, and limitations on Anopheles vector management and surveillance 29 . In Africa, the patterns of malaria incidence are widespread especially in sub-Saharan desert areas by human migration and geographical changes such as dam construction, while malaria transmission in ROK is confined specifically to northern border areas 30 . This suggests that infected mosquitoes migrating from North Korea may be one of the major factors of malaria incidence in ROK, considering the restriction of human migration in the northern border areas of ROK. Therefore, it is important to trace the originality of malaria mosquitoes and their microbiome in hyperendemic areas by employing studies such as aerial netting methods to collect migrating mosquitoes at high altitudes from the north of DMZ, which could help to eradicate malaria in ROK 31 . The microbiome of the hyperendemic and hypoendemic regions has high dissimilarity (Fig. 6). Nevertheless, this study faces certain limitations in comprehensively understanding the regional variation of mosquito microbiome within the hyperendemic area. Even though the hyperendemic 2-M sample was collected from a hyperendemic area, it had a microbiome that was similar to that found in samples collected from a hypoendemic area. Moreover, samples from the hyperendemic area showed three clustering, although this did not seem to be strongly related to the distance between each region. Further details on the environment of the region where the mosquito was collected (such as the mosquito's aquatic habitat, blood sources, etc.) as well as whether the mosquito is Plasmodium-infected will be needed to comprehend the variation by region.
Interestingly, our study demonstrates that Pseudomonas is significantly different between the mosquito microbiome profiles of hyperendemic and hypoendemic areas, showing Pseudomonas was dominant in hypoendemic areas. Pseudomonas is commonly found in the gut microbiota of Anopheles mosquitoes [32][33][34][35] , indicating that Pseudomonas is capable of efficiently adapting to the mid-gut environment of Anopheles mosquitoes. In addition, a previous study has pinpointed that Pseudomonas is the most prevalent microbiota in the Plasmodium-negative groups 36 and protects mosquitoes from the invasion of malaria parasites 37 . Interestingly, our study shows a low proportion of the Psuedomonas population of microbiome profiles in the hyperendemic areas, indicating that www.nature.com/scientificreports/ there might be some factors such as malaria parasites to disturb the balance of microbiota. It has been reported that Pseudomonas synxantha inhibits the growth of Mycobacterium tuberculosis 38 causing tuberculosis, while little is known about the role of Pseudomonas synxantha in the mosquito. Having said that, it will be intriguing to elucidate the function of Pseudomonas which is the most dominant microbiome in hypoendemic areas in our study. Another focus of our study was the identification of species biomarkers specific to different endemic areas of malaria incidence in ROK. Therefore, further studies on spatial-temporal patterns of the mosquito microbiome Anopheles sinensis microbiome diversity according to epidemiology and organ This figure offers images that were shot in a multiple ranges of angles. The PCoA plot is a graph in which the distance between areas is estimated based on their dissimilarity in terms of relative abundance (Bray-Curtis). It can be seen that samples from the hyperendemic area are spread by a greater distance than those from the hypoendemic area, while samples from the hypoendemic area are closer together. www.nature.com/scientificreports/ are likely to be important to determine a seasonal variation of microbiome aligned with malaria incidence as well as tracing mosquito migration from elsewhere. Although the major function of the mid-gut and salivary glands of mosquitoes is to digest diverse nutrients, these organs are crucial for the growth of gut-associated microbiota 39,40 and pathogen transmission of mosquito-borne diseases [41][42][43][44] . The transmission of malaria parasites is dependent on maintaining parasite growth in the mid-gut and culminating virulent sporozoite stages in the salivary glands 45 . Therefore, the region-specific microbiome such as Staphylococcus, Erwinia, Enterobacteriaceae, Serratia, Pantoea, and Acinetobacter identified in each part of the hyperendemic area of this study has a potential role to interact with malaria parasites. Serratia marcescens is known to block the sporogonic development of P. vivax parasites in An. albimanus 46 , and Acinetobacter species increase the resistance of An. gambiae to Plasmodium development partly by the induction of anti-Plasmodium factors in Imd pathway 47 . Therefore, the potential roles of these region-specific microbiotas in the hyperendemic areas on the effect of malaria transmission remain to be investigated. In addition to the effects of the microbiota on parasite development, the microbiota is known to modulate host behaviors and thus increase vector competence. It has been reported that microbiome-Gut-Brain-Axis communication can change the sensory perception ability and blood-sucking behaviors of mosquitoes by specific neuropeptides secreted from certain microbiota 48 . Mosquitoes infected with malaria parasites show stronger blood-sucking behaviors, increasing the number of biting and precise finding of human hosts 49 , indicating that there might be some changes in microbiome-Gut-Brain-Axis communication in modified microbiota community. Interestingly, it has been demonstrated that the rapid proliferation of Pseudomonas in the mosquito mid-gut after blood-feeding induces the secretion of serotonin to suppress appetite, followed by the secretion of neuropeptide Y (NPY) in the mosquito brain 48,50,51 , which has pivotal roles in modulating sensory sensitivity such as olfaction 51 . Previous data indicate that one of the neuropeptides, tachykinin, also modulates olfactory pathways and insulin pathways in insects 52 and synaptic contacts between terminals of tachykinergic neurons and NPY neurons regulate NPY neuronal activity in rats 53 , suggesting that neuropeptides secreted from microbiome can alter host finding behaviors of the mosquito by these neuronal pathways. Taken together, our study provides evidence that the microbiome composition in the mid-gut and salivary glands show spatially dynamic patterns, and it may influence the malaria incidence in ROK. Additional spatio-temporal variation of the microbiome will be providing more precise clues to understanding the origin of malaria mosquitoes in hyperendemic areas. It will be also interesting to elucidate the molecular and neuronal mechanisms of microbiome-Gut-Brain-Axis communication among host-microbiome-parasite underlying the modulation of host-finding behaviors. Although this study provides a snapshot of microbiome patterns in malaria-endemic areas in ROK, future studies will be required to estimate the Plasmodium concentration inside the mosquitoes samples to fully understand the exact patterns and dynamics of the microbiome in an individual mosquito, which could explain the correlation between Plasmodium spores and microbiome composition. Later on, fully understanding of microbiome patterns of malaria mosquitoes in the ROK and Korean peninsula will mitigate the efforts of the malaria eradication program.

Methods
Sample collection and identification. Female adult Anopheles sinensis mosquitoes were collected in 2020 (June 22 to 26) from 12 rural regions in ROK with the assistance of the regional center for vector surveillance against climate change supported by the Korea Centers for Disease Control and Prevention 54 . The collection sites were categorized as hyperendemic area if the number of malaria cases per 100,000 persons exceeds 1 patient. Otherwise, it was categorized as hypoendemic area (Fig. 1). The study regions were determined using statistics of malaria patient. The analysis could have improperly accounted for habitat or ecology since each region was not ecologically analyzed. Mosquitoes were sorted into genus using the morphological keys 55 , and species identifications were confirmed using a diagnostic PCR assay based on DNA barcode analysis 56 . Mosquito body parts such as legs and wings were transferred to a 1.5ml tube (Axygen, USA), after which DNA extraction was carried from mosquito parts using G-spin™ Total DNA Extraction kit (iNtRON's, Korea) according to the manufacturer's instructions. PCR cycle parameters involved an initial denaturation at 95 °C for 3 min, 35 cycles of 30 s at 95 °C, 30 s at 63 °C, and 2 min at 72 °C. A final extension at 72 °C for 10 min was completed. PCR products were subjected to electrophoresis on a 1.5% agarose gel and visualized under ultraviolet light. DNA extraction. The DNA extraction targeted Adult female mosquito mid-gut and salivary gland. Before to dissection, mosquitoes were surface sterilized for 1 min in 70% ethanol, then dissected in PBS (Phosphate-Buffered Saline). The dissecting stereo microscope working area was likewise sanitized with 70% ethanol during the dissection. mid-guts and salivary glands were pooled (each sample comprised of 5 mosquitoes) and kept at − 80 °C. Whole genome DNA was extracted under aseptic conditions using the Power soil kit (Qiagen, Germany) according to the manufacturer's recommendations. Samples were used for 16S metagenomics analysis after DNA quality and quantity were checked. Two 16S rRNA-tagged libraries based on amplicons were created. www.nature.com/scientificreports/ paired-end reads. The paired-end sequencing reads were first imported into the QIIME2 v.2021.11 57 pipelines for 'quantitative insights into microbial ecology' , and the pipeline's v.2021.11 was used to process and analyze the sequencing reads. The denoise-paired command was used to rectify errors, eliminate chimeras, and integrate paired-end reads using the DADA2 plugin in QIIME 2 58 . The sequence generated was utilized to compare the bacterial composition and taxonomy analysis later on. Moreover, chunlab analytical pipeline PKSSU 4.0 DB was used to group sequences at a 3% divergence (or 97% similarity) level, sequences were denoised and allocated to an operational taxonomic unit (OTU). Chimeric sequences that could not be definitively linked to an OTU were eliminated. OTUs that could not be taxonomically categorized were labeled "unclassified" and were excluded from further analysis. Phylum, genus, and species level data-sets were utilized for our analysis. At each level, count and ratio data are provided equally, and ratio data was used for analysis. For comparison, the ratio data were subdivided into epidemiology and organ ( Supplementary Informations 2-4).

16S sequencing and taxonomic analysis.
Bacterial profiles. At the phylum, genus, and species levels, read count and abundance data for bacterial OTUs were evaluated. Low abundance taxa with a value of less than 1% were grouped into an "Other bacteria" category. EzBioCloud 16S database is designed to be best performed for species-level identification even though there is a limitation due to the lack of sequence differences in some closely related species. The combination of EzBioCloud 16S DB and sensitive bioinformatics pipelines allows us a species-level exploration of microbiome data.
Species biomarker. Species biomarkers were used to find microbiota specific to each collection site. All analyses were conducted using "16S-based Microbiome Taxonomic Profiles (MTP)-Comparative Analyzer for MTP sets-Taxonomic biomarker discovery", using the EzBioCloud application. A group containing samples from all regions and a group containing only samples from a specific region were compared using the Kruskal-Wallis H test, and among the microbiota with statistically significant values, the two or three microbiotas with the highest ratio were selected as species biomarkers ( Supplementary Information 6).
Statistical analysis. The Shannon diversity, Jaccard, and the Bray-Curtis dissimilarity index were used to analyze microbial diversity within (alpha diversity) and between (beta diversity) samples in QIIME2. The average Shannon indices that resulted were reported and compared between samples using pairwise Kruskal-Wallis tests with Benjamini-Hochberg FDR corrections for multiple comparisons. The results of the Jaccard and Bray-Curtis dissimilarity indices were compared between samples using paired PERMANOVA tests (999 permutations) with FDR corrections. The statistical analyses described above were performed using EZbiocloud 59 (https:// www. ezbio cloud. net/) and QIIME2 60 .

Data availability
The metagenome reads obtained from this study have been deposited in NCBI under the BioProject PRJNA844511. All relevant data are either within the paper or will be submitted and will be available in a public repository at NCBI.